*** the approach here is:
* (1) within each combination of nuisance-factor and stratum, any observation (regardless of whether they chose cap, no cap, or did not participate) could have been randomized into treatment with equal probability --> see RANDOMNESS 1 part
* (2) even conditional on the nuisance-factor assignment, difference nuisance-factor/stratum groups could have received somewhat varying numbers of treatment, albeit not independently because the randomization was blocked and stratified --> see RANDOMNESS 2 part
* NB: this is implicitly conditioning on the nuisance factor assignment not only of the judges who participated but also those who didn't. The assignment mechanism could probably have created a few other combinations of treatment assignment compatible with the participating judges' nuisance factor assignments but not the non-participating ones.

tempfile alldata data crand_2018_a_location crand_2018_sympathy rerandomized
cap log close
log using analysis_output/randomization_tests, replace

*** things to loop over
foreach treatment_factor in "sympathy" "a_location" {
	if "`treatment_factor'" == "a_location" local nuisance2 "sympathy"
	else local nuisance2 "a_location"
	local nuisance_factors "forum `nuisance2'"
	local factors "`nuisance_factors' `treatment_factor'"
	di "TREATMENT: `treatment_factor'"
foreach uncertainjudge in "District" "Magistrate" { // this is for 2017 randomization where one WNP judge could be either District or Magistrate -- matters only for location
	di "uncertain judge --> `uncertainjudge'"
foreach exclusion in "" `"if Q4!="No""' "if prior!=1" `"if Q4!="No" & prior!=1"' {
	di `"exclusion: `exclusion'"'
	use alldata.dta `exclusion', clear
	qui save `alldata', replace

qui {
**************************************************************************************************************
******************************************       2018      ***************************************************
**************************************************************************************************************

use `alldata' if year==2018, clear
assert inlist(cap,0,1) // confirming 0/1 coding of cap, which is essential for everything below
collapse (sum) k=cap (count) n=cap, by(`nuisance_factors' judge)
save `data', replace

use raw_data/2018/randomization_results_as_used_2018_anonymous.dta, clear
do labels
encode accident_state, gen(a_location) label(location)
rename (forum sympathy judgetype) =_
encode forum, g(forum) label(forum)
encode sympathy, g(sympathy) label(sympathy)
encode judgetype, g(judgetype) label(judge)
recode `treatment_factor' (-1=0)
collapse (count) N = `treatment_factor' (sum) T = `treatment_factor', by(`nuisance_factors' judgetype)
merge 1:1 `nuisance_factors' judge using `data', assert(1 3) keep(3) nogenerate // esp. Circuit judges were few so 3 out of 4 cells for them are empty
compress
save `data', replace

levelsof judge, local(judgetypes)


*** RANDOMNESS 1: for given assignment vector, randomly assigning them to judges who (i) did not participate (-->f), (iia) applied the cap or (iib) not (-->g)

local z = 0 // for labelling the nuisance-specific datasets
forvalues forum=-1(2)1 {
forvalues nuisance=-1(2)1 {
	tempfile data`++z'
	local otherdata = 0
	foreach judge of local judgetypes {
		use `data' if judge==`judge' & `nuisance2'==`nuisance' & forum==`forum', clear
		if _N==0 continue // if stratum is empty (i.e., no judge of this type participated for this forum-nuisance combination), move on to next stratum
		else {
			qui {
				* draw at random t assigned treatments for n participating judges from treatment assignments for N registered judges, T of which were treated) ...
				gen min = max(0,T-N+n)
				gen max = min(n,T)
				expand max-min+1
				gen byte t = min + _n-1
				gen double f = hypergeometricp(N,T,n,t)
				recode f (.=1) // Stata generates missing values when there is only one possible allocation of successes
				
				* ... times, within that, distribution of k_t caps when drawing t outcomes fron n participating judges, k of which applied cap
				replace min = max(0,k-n+t)
				replace max = min(t,k)
				expand max-min+1
				bys t: gen byte k_t = min + _n-1
				gen double g = hypergeometricp(n,k,t,k_t)
				recode g (.=1)
				
				* --> joint distribution
				gen double p = f*g
				keep p `nuisance_factors' n k t k_t
			}
			
			if `otherdata++'>0 { // combine with other strata, if any -- can ++ here because we only get here if "continue" above wasn't triggered, i.e., we have data here
				rename * =_2
				cross using `data`z'' // this is the data from the other strata processed so far in this nuisance cell
				replace p = p * p_2
				foreach var of varlist n k t k_t {
					replace `var' = `var' + `var'_2
				}
				collapse (sum) p, by(`nuisance_factors' n k t k_t)
			}
			
			save `data`z'', replace
		}
		}
	
	use `data`z'', clear
	rename * =_`z'
	save `data`z'', replace
	
}
}


*** RANDOMNESS 2: assignment vector could have been different, albeit only in limited ways when conditioning on nuisance factors
	* balancing of factor level assignment means flipping treatment for one nuisance factor combo necessarily flips observations
		* (a) with the same nuisance factors (i.e., all obs in the same z group) AND
		* (b) the respective opposite that balances the rows (that's z groups 1<-->4 and 2<-->3)
		** Example: flipping 000 to 001 flips (a) 001 to 000 AND (b) 111 to 110 and 110 to 111
		* ---> so these can only be flipped together!
	* flipping is a mere relabeling of treatment and control, so we get the resulting distribution by changing t to n-t and k_t to k-k_t
	* if 2018 were only data, we would need to enumerate only one half of the possible flip combos because the other is the same except treatment multiplied by -1 --> same |t|
		* however, that doesn't work when combining with 2017 data

tempfile data14
foreach pair in 14 23 { // those are the nuisance factor combinations that can only be flipped together (e.g., 1 and 4)
	local a = floor(`pair'/10)
	local b = mod(`pair',10)
	use `data`a'', clear
	cross using `data`b''
	gen p`pair' = p_`a' * p_`b' / 2 // division by two anticipates doubling of data below
	drop p_*
	expand 2, gen(copy)
	foreach x in `a' `b' {
		replace t_`x' = n_`x' - t_`x' if copy // this and next line "mirror" the data -- since treatment and control could have been reversed
		replace k_t_`x' = k_`x' - k_t_`x' if copy
		}
	drop copy
	qui ds p`pair', not
	collapse (sum) p`pair', by(`r(varlist)')
	if `pair'==14 save `data14', replace
}

cross using `data14'
gen p = p14*p23
drop p14 p23

save `crand_2018_`treatment_factor'', replace

**************************************************************************************************************
***************************************    Sympathy 2017 ******************************************************
**************************************************************************************************************

if "`treatment_factor'"=="sympathy" { // sympathy randomization was last and not stratified in 2017 -- so we also don't need to worry about missing judgetype (3 judges)!
	
	use `factors' cap year using `alldata' if year==2017, clear
	collapse (sum) k=cap (count) n=cap, by(`nuisance_factors')
	save `data', replace

	use raw_data/2017/randomization_result_anonymous_as_used_on_20170410.dta, clear
	do labels
	encode w_forum, g(forum) label(forum)
	encode w_accidentstate, g(a_location) label(location)
	encode w_sympathy, g(sympathy) label(sympathy)
	recode `treatment_factor' (-1=0)
	collapse (count) N = `treatment_factor' (sum) T = `treatment_factor', by(`nuisance_factors')
	assert a_location=="K":location & forum=="S":forum if mod(N,2)!=0 // only nuisance factor combo with uneven N is this one -- this matters for randomness 2 below
	merge 1:1 `nuisance_factors' using `data', assert(1 3) keep(3) nogenerate
	compress
	save `data', replace
	
	local z = 0 // for labelling the nuisance-specific datasets
	forvalues forum=-1(2)1 { // this and next loop need to be set up identically to 2018 loops to match the nuisance cells
	forvalues nuisance=-1(2)1 {
		tempfile data`++z'
		use `data' if `nuisance2'==`nuisance' & forum==`forum', clear
		
		* this part is the same as Randomness 1 in 2018 but without the complications of judge stratification (which wasn't done for SYMPATHY in 2017)
		qui {
			* draw at random t assigned treatments for n participating judges from treatment assignments for N registered judges, T of which were treated) ...
			gen min = max(0,T-N+n)
			gen max = min(n,T)
			expand max-min+1
			gen byte t = min + _n-1
			gen double f = hypergeometricp(N,T,n,t)
			recode f (.=1) // Stata generates missing values when there is only one possible allocation of successes
			
			* ... times, within that, distribution of k_t caps when drawing t outcomes fron n participating judges, k of which applied cap
			replace min = max(0,k-n+t)
			replace max = min(t,k)
			expand max-min+1
			bys t: gen byte k_t = min + _n-1
			gen double g = hypergeometricp(n,k,t,k_t)
			recode g (.=1)
			
			* --> joint distribution
			gen double p = f*g
			keep p `nuisance_factors' n k t k_t
		}
		
		* this part is the equivalent of Randomness 2 in 2018: if N is uneven, T could have been (N+1)/2 or (N-1)/2 -- equivalently, switch T to N-T, t to n-t, and k_t to k-k_t
		* this only concerns forum=S X a_loc=K" combo since all others have even N (see assert above when assembling data)
		* hence we don't need to account for the dependence across nuisance cells from balancing aspect of code (which tried to keep imbalance to at most 1)
		
		if `forum'=="S":forum & `nuisance'=="K":location {
			expand 2, gen(copy)
			replace t = n - t if copy // this and next line "mirror" the data -- since treatment and control could have been reversed
			replace k_t = k - k_t if copy
			drop copy
			qui ds p, not
			collapse (sum) p, by(`r(varlist)')
			replace p = p/2
			}
		
		* cleanup
		rename * =_`z'
		save `data`z'', replace	
		
	}
	}

	use `data1', clear
	rename p_1 p
	forvalues x=2/4 {
		cross using `data`x''
		replace p = p*p_`x'
		drop p_`x'
	}
}


**************************************************************************************************************
***************************************   Location 2017 ******************************************************
**************************************************************************************************************

if "`treatment_factor'"=="a_location" { 
	
	*** Forum = WY
		** method notes: there are so few possibilities with the 5 obs in 2017 that it seems easier, less error-prone, and more transparent to enumerate them explicitly
	use `factors' judge cap year using `alldata' if year==2017, clear
	keep if forum=="W":forum
	
	assert _N==5
	assert judge!="Circuit":judge // no Circuit judges
	assert cap == 1 if judge=="Bankruptcy":judge // Bankruptcy judges are there but decided the same, so treatment assignment irrelevant	
	replace judge="`uncertainjudge'":judge if mi(judge) // the one participant with missing judge type who could be either district or magistrate (see data assembly code)
	
	tempfile W_2017_location
	foreach district_switch in 1 -1 { // 1 changes nothing, -1 reverses treatment assignment (which was equally likely) -- no other assignments were possible because of balanced assignment within judge type
	foreach magistrate_switch in 1 -1 { // id.
		preserve
		replace a_location = a_location * `district_switch' if judge=="District":judge
		replace a_location = a_location * `magistrate_switch' if judge=="Magistrate":judge
		recode a_location cap (-1 = 0) // cap is probably coded (0,1) so far but just in case
		gen k_t_ = cap*a_location
		collapse (count) n_ = cap (sum) t_ = a_location (sum) k_ = cap (sum) k_t_, by(`nuisance_factors')
		assert forum==-1 & inlist(sympathy,-1,1)
		gen z = 1 if sympathy==-1 // cf. z labelling in 2018 Randomness 1
		replace z = 2 if sympathy==1
		rename (forum sympathy) =_
		gen helper = 1
		reshape wide n_ t_ k_ k_t_ sympathy forum, i(helper) j(z)
		drop helper
		gen double p = 1/4
		cap append using `W_2017_location'
		save `W_2017_location', replace
		restore
	}
	}
	
	
	*** Forum = SD
		** method notes:
			* STEP 1: get sympathy-conditional distribution of treatments across judge-sympathy cells (! note the sympathy conditioning: without it, it would just be half of each judge cell)
			* by Bayes' Rule, pr(t|S)=pr(S|t)pr(t)/pr(S). We can calculate pr(S|t) and pr(t) straightforwardly from the assignment mechanism. If we do this for all t that can yield S, we can sum them to obtain pr(S)
			* STEP 2: assign these treatments randomly to participants and non-participants within these judge-sympathy cells
			* this step is similar to before: hypergeometricp(# judges in randomization, # assigned treatment, n=# who participated, n_t=# treated & participated)*hypergeometricp(n,k=caps,n_t,k_t=caps among treated)
	
		* STEP 1
			* get raw inputs, transform into frequencies: # randomized judges by type and sympathy assignments (because we need to condition on those)
			use raw_data/2017/randomization_result_anonymous_as_used_on_20170410.dta if w_forum=="S", clear
			do labels
			encode judgetype, g(judge) label(judge)
			levelsof judge, local(judgetypes)
			encode w_sympathy, g(sympathy) label(sympathy)
			assert "S":sympathy==1 // basis for coding etc. later
			recode sympathy (-1 = 0) // required for suffixing later
			count if sympathy==1
			scalar S = r(N)
			scalar N = _N
			contract sympathy judge, freq(n_)
			
			reshape wide n_, i(sympathy) j(judge)
			rename n_? =_
			gen helper = 1 // only used to facilitate the next line: reshape, which requires argument for i()
			reshape wide n_?_, i(helper) j(sympathy)
			drop helper
			foreach judge of local judgetypes {
				egen n_`judge' = rowtotal(n_`judge'_?)
			}
			
			* A1: possible distribution of # treated (nt) within and across judge types
			foreach judge of local judgetypes {
				gen nt_`judge' = n_`judge'/2 // treatment balance within judge cells means nt = n/2 ...
				expand 2 if mod(nt_`judge',1)!=0, gen(expander) // unless n is uneven, in which case +/- 1/2
				replace nt_`judge' = nt_`judge' + (2*expander-1)/2 if mod(nt_`judge',1)!=0 // "if" should be unnecessary but just in case
				drop expander
			}
			egen nt = rowtotal(nt_*)
			drop if !inrange(nt,N/2-1/2,N/2+1/2) // the randomization code enforced a global balance such that total imbalance of treated vs. not treated could be a most +/- 1
			scalar G=_N // this scalar keeps track of how many big forks we considered, which were all equally likely -- later we divided p by this scalar
			assert G==6 // that's because we had three groups with uneven participant numbers -- but if we ran this on different data, this line should be dropped
			compress
			
			* A2: possible values of # treated sympathetic for all judge types combined (nt_s) 
			gen int nt_s = floor(nt/2) // half of each group (treated, untreated) is assigned to sympathetic (S), and if it's the only uneven group, the one leftover is assigned N. But ...
			if mod(nt,2)!=0 & mod(N-nt,2)!=0  { // ... if both groups have a leftover, then they are randomly assigned to N and S
				expand 2, gen(plus)
				replace nt_s = nt_s + plus // plus should evaluate for 0 to the original and 1 for the duplicated observation
				drop plus
				scalar G = G*2 // imbalance was equally likely to have one more sympathetic or unsympathetic treated
			}
			
			* A3: possible division of nt_s between judge types
			foreach judge of local judgetypes {
				gen id = _n
				expand min(nt_`judge', n_`judge'_1) + 1 // +1 because nt_?_1 may be zero (even if nt=1)
				bys id: gen nt_`judge'_1 = _n - 1  // all possibilities of # judges of this type who have sympathy=1 and are treated
				gen nt_`judge'_0 = nt_`judge' - nt_`judge'_1 // this is for later when we assign t to cap, no cap, or missing
				drop if nt_`judge'_0 > n_`judge'_0
				drop id
			}
			egen nt1 = rowtotal(nt_?_1)
			drop if nt1!=nt_s
			compress
			
			* P1*P2: unconditional probability of observing this distribution of treatment assignments across judge types
			gen double p = 1/G // see above: distribution of extra treatment and extra sympathy across judge types	
			
			* P3: probability of observing this distribution of S conditional on treatment assignments across judge types
			gen nnt = N - nt // need to keep track of non-treated units, and probability of their S assignment -- we are conditioning on all S assignments!
			gen nnt_s = S - nt_s
			foreach judge of local judgetypes {
				replace p = p*min(1,hypergeometricp(nt,nt_`judge',nt_s,nt_`judge'_1)) * /// min because hypergeometricp evaluates to . if p=1
							min(1,hypergeometricp(nnt,n_`judge'-nt_`judge',nnt_s,n_`judge'_1-nt_`judge'_1))
				replace nt = nt - nt_`judge'
				replace nt_s = nt_s - nt_`judge'_1
				replace nnt = nnt - (n_`judge'-nt_`judge')
				replace nnt_s = nnt_s - (n_`judge'_1-nt_`judge'_1)
			}
			
			sum p
			replace p = p/r(sum) // the denominator of Bayes' formula
			
			keep n_?_? nt_?_? p
			compress
			tempfile 2017_st
			save `2017_st', replace
		
			
		* STEP 2
			* get inputs: actual participant numbers and # cap (by sympathy and judgetype)
			use `alldata', clear
			keep if forum=="S":forum & year==2017 & !mi(judgetype)
			recode sympathy (-1 = 0)
			forvalues s=0/1 {
				count if sympathy==`s'
				scalar n_`s'=r(N)
				count if sympathy==`s' & cap==1
				scalar k_`s'=r(N)
			}
			collapse (count) np_ = cap (sum) k_ = cap, by(judge sympathy)
			reshape wide np_ k_, i(sympathy) j(judge)
			rename (np_? k_?) =_
			gen helper = 1
			reshape wide n* k*, i(helper) j(sympathy)
			drop helper
			assert _N==1
			foreach j of local judgetypes {
			forvalues s=0/1 {
				scalar np_`j'_`s' = np_`j'_`s'[1]
				scalar k_`j'_`s' = k_`j'_`s'[1]
			}
			}
			
			* calculate joint probabilities of assignments of treatments to participants (who did or did not apply cap)
			use `2017_st', clear // need to keep the lines of these data together because the randomization across judges, sympathy isn't independent
			foreach j of local judgetypes { // cf. code for Randomness 1 in 2018
			forvalues s=0/1 { 
				gen i = _n	
				expand min(np_`j'_`s',nt_`j'_`s')-max(0,nt_`j'_`s'-n_`j'_`s'+np_`j'_`s')+1 // t among np: <=min(np,t), >=max(0,t-attrition) (a=n-np)
				bys i: gen byte t_`j'_`s' = max(0,nt_`j'_`s'-n_`j'_`s'+np_`j'_`s') + _n - 1 // t (nt) is # treated among participated (randomized)
				gen double f = min(1,hypergeometricp(n_`j'_`s',nt_`j'_`s',np_`j'_`s',t_`j'_`s'))
				replace i=_n
				expand min(t_`j'_`s',k_`j'_`s')-max(0,k_`j'_`s'-np_`j'_`s'+t_`j'_`s')+1 // k among t: <=min(t,k), >=max(0,k-np+t)
				bys i: gen k_t_`j'_`s' = max(0,k_`j'_`s'-np_`j'_`s'+t_`j'_`s') + _n - 1
				gen double g = min(1,hypergeometricp(np_`j'_`s',t_`j'_`s',k_`j'_`s',k_t_`j'_`s'))
				replace p = p*f*g
				sum p, meanonly
				assert float(r(sum))==1
				drop f g i
			}
			}
			
			* collapse over judgetype (no longer needed)
			forvalues s=0/1 {
				egen t_`s'	 = rowtotal(t_?_`s')
				egen k_t_`s' = rowtotal(k_t_?_`s')
				drop t_?_`s' k_t_?_`s'
			}
			collapse (sum) p, by(t_? k_t_?)
			compress
			
	*** joining the two fora, and then 2018
	
		* clean-up of S in 2017: add common information including identifiers
		forvalues s=0/1 {
			local z = `s' + 3 // cf. z labelling in 2018 Randomness 1
			rename t_`s' t_`z'
			rename k_t_`s' k_t_`z'
			gen n_`z' = n_`s'
			gen k_`z' = k_`s'
			gen forum_`z' = 1
			gen sympathy_`z' = 2*`s'-1
		}
		compress
		
		* join S and W in 2017
		rename p p_S
		cross using `W_2017_location'
		replace p = p*p_S
		drop p_S
		
	}
	
	
**************************************************************************************************************
****************************   cross with 2018   **************************************************************
**************************************************************************************************************

rename * =_2017
cross using `crand_2018_`treatment_factor''

* check that group labelling is consistent
forvalues z=1/4 {
	assert forum_`z'==forum_`z'_2017
	assert `nuisance2'_`z'==`nuisance2'_`z'_2017
}

* combine 2017 and 2018 values
replace p = p*p_2017
forvalues z=1/4 {
	foreach stat in n k t k_t {
		replace `stat'_`z' = `stat'_`z' + `stat'_`z'_2017
	}
}

drop *_2017
qui ds p, not
collapse (sum) p, by(`r(varlist)')

* test block
sum p, meanonly
assert float(r(sum)) == 1
forvalues x=1/4 {
	assert n_`x'>=k_`x'
	assert n_`x'>=t_`x'
	assert k_t_`x'<=min(k_`x',t_`x')
}	

save `rerandomized', replace
			

**************************************************************************************************************
******************************   test     ********************************************************************
**************************************************************************************************************

use `alldata', clear
if "`treatment_factor'" == "a_location" { // for sympathy, judge type doesn't matter, i.e., all obs. are useable -- see note at top of 2017 sympathy randomization
	local condition `"if mi(judgetype) & year==2017 & forum=="W":forum & sympathy=="N":sympathy"'
	qui count `condition'
	assert r(N)==1
	replace judgetype="`uncertainjudge'":judge `condition'
	drop if mi(judge)
	}

* verify numbering (--> from 2018, Randomness 1)
preserve
	use `rerandomized', clear
	run labels.do
	assert "W":forum==-1
	assert forum_1 == forum_2
	assert forum_3 == forum_4
	assert forum_1 =="W":forum
	assert forum_3 =="S":forum
	assert `nuisance2'_1 == `nuisance2'_3
	assert `nuisance2'_2 == `nuisance2'_4
	assert `nuisance2'_1 == -1
	assert `nuisance2'_2 == 1
restore

* generate count statistics (nomenclature: upper case for observed data, lower case for randomization distribution)
local z = 0 
forvalues forum=-1(2)1 { // ordering as verified just now
forvalues nuisance=-1(2)1 {
	local `z++'
	count if `treatment_factor'==1 & forum==`forum' & `nuisance2'==`nuisance'
	scalar T_`z' = r(N)
	count if `treatment_factor'==1 & forum==`forum' & `nuisance2'==`nuisance' & cap==1
	scalar K_t_`z' = r(N)
}
}
scalar T_W = T_1 + T_2
scalar K_t_W = K_t_1 + K_t_2
scalar T_S = T_3 + T_4
scalar K_t_S = K_t_3 + K_t_4

* generate difference in proportion
foreach forum in W S {
	prtest cap if forum=="`forum'":forum, by(`treatment_factor')
	scalar D_`forum' = -r(P_diff) // "-" because prtest subtracts the second group (1 -- sympathetic defendant / Nebraska) from the first (-1 -- unsympa / K)
}
if "`treatment_factor'" == "a_location" scalar D_W = -D_W // getting sign right matters for one-sided tests and, for location, comparisons between fora

* load randomization distribution and generate its count and proportion statistics for comparison
use `rerandomized', clear

gen k_t_W = k_t_1+k_t_2
gen t_W = t_1 + t_2
gen k_t_S = k_t_3+k_t_4
gen t_S = t_3 + t_4

gen d_W = (k_t_W/t_W)-(k_1+k_2-k_t_W)/(n_1+n_2-t_W)
gen d_S = (k_t_S/t_S)-(k_3+k_4-k_t_S)/(n_3+n_4-t_S)
if "`treatment_factor'" == "a_location" replace d_W = -d_W // cf. parallel operation for actual data a couple lines above

}

* within fora comparison of count stats
foreach forum in W S {
	di "`forum': " K_t_`forum' " successes, diff in proportions = " %5.4f D_`forum'
foreach restriction in  "forum-wide" "cell-by-cell" {
	preserve
	
	* condition on # treated units (if not, code: "if "`forum'" == "W" keep p *_1 *_2 *_W; else keep p *_3 *_4 *_S")
	if "`restriction'"=="forum-wide" qui keep if t_`forum'==T_`forum'
	if "`restriction'"=="cell-by-cell" & "`forum'"=="W" qui keep if t_1==T_1 & t_2==T_2
	if "`restriction'"=="cell-by-cell" & "`forum'"=="S" qui keep if t_3==T_3 & t_4==T_4
	sum p, meanonly
	qui replace p = p/r(sum) // necessary with conditioning because we have excluded many possibilities; without restrictions it's redundant
	
	* sum probabilities of more "extreme" events
	collapse (sum) p, by(k_t_`forum' d_`forum') // two ways: number of "successes" (k) or diff. in proportion (d) (identical for simple Fisher test)
	isid k_t_`forum' // logically has to be the same
	sum p if k_t_`forum' == K_t_`forum', meanonly // this only serves to get the cutoff probability
	assert r(min)==r(max) // see previous line -- should be a single value
	sum p if float(p)<=float(r(min)), meanonly // sum of probabilities of all less likely events
	scalar p_prob = r(sum)
	sum p if float(abs(d_`forum')) >= float(abs(D_`forum')), meanonly
	assert p_prob == r(sum), rc0 // p-value should be same whether "extreme" defined by prob. or diff. in proportion (algebraically guaranteed for simple Fisher test)
	di "p = " %5.4f r(sum)  _continue
	sum p if float(d_`forum') >= float(D_`forum'), meanonly
	di " ("  %5.4f r(sum) " one-sided) (`restriction')"
	restore
}
}

* between fora comparison of difference in proportions
di "Between fora: difference in proportion"
foreach restriction in "" "forum-wide" "cell-by-cell" {
	if "`restriction'"!="" di "restriction to equal number of treated units: `restriction'"
foreach weight in "unweighted" "variance-weighted" {
	preserve
	if "`restriction'"=="forum-wide" qui keep if t_W==T_W & t_S==T_S
	if "`restriction'"=="cell-by-cell" qui keep if t_1==T_1 & t_2==T_2 & t_3==T_3 & t_4==T_4
	sum p, meanonly
	qui replace p = p/r(sum) // necessary with restrictions because we have excluded many possibilities; without restrictions it's redundant
	if "`weight'"=="variance-weighted" {
		qui sum d_W // [iw=p] generates error because the sum command calculates a negative variance ...
		scalar V_W = r(Var)
		qui sum d_S // [iw=p]
		scalar V_S = r(Var)
		scalar l = V_S/(V_W+V_S)
		qui gen d = l*d_W - (1-l)*d_S
		scalar D = l*D_W - (1-l)*D_S
	}
	else {
		qui gen d = d_W - d_S
		scalar D = D_W - D_S
	}
	sum p if float(abs(d))>=float(abs(D)), meanonly
	di "   `weight' diff in proportions = " %4.3f D ", p = " %5.4f r(sum) _continue
	sum p if float(d)>=float(D), meanonly
	di " ("  %5.4f r(sum) " one-sided)"
	restore
}
}

}
}
}

log close